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ABSTRACT 

Simulation of fluid flow in porous media using lattice Boltzmann method depends 
on how effectively collision, streaming and boundary conditions are implemented at 
micro level in view of the macroscopic behaviour of fluid. While collision and 
streaming have been extensively researched and attained effective formulation, 
boundary conditions still remains to be studied thoroughly. This paper studies various 
boundary conditions that are defined to simulate non-Newtonian fluid flow in porous 
media using lattice Boltzmann method. Further, these conditions are applied to 
simulate the problem of non-Newtonian forced convection in porous media and the 
variation inflow regimes and rate of heat transfer is studied based on the variation in 
velocity and thermal boundary behavior. Though velocity boundary conditions did not 
produce any difference in the flow regimes, thermal boundary conditions produced 
significant variation in rate of heat transfer. 
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1. INTRODUCTION 

Implementation of boundary conditions is one of the most important steps in simulating fluid 
flow, which brings stability and accuracy to the numerical method being used. Because of the 
kinetic approach of Lattice Boltzmann method (LBM) and macroscopic behavior of fluid flow 
not being very sensitive to microscopic particle dynamics, implementing boundary conditions 
becomes easier as compared to macroscopic methods. It is this aspect of LBM which makes it 
a preferred numerical tool over others to simulate fluid flow through complex geometries, 
especially when there is a particle-wall interaction like in the case of porous media, colloids 
and confined flows. At the same time, kinetic approach deals with fluid behavior at a much 
smaller scale compared to macroscopic scale, therefore, fluid wall interaction has to be 
incorporated carefully. Generally, the velocity space in a kinetic approach to boundary 
conditions is split two parts: incoming velocities and outgoing velocities (Gross et al. (1957), 
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Gross and Jackson (1959)), which are distinguished based on discontinuity in the distribution 
function. Different approaches to consider boundary conditions have been proposed (Walsh et 
al. (2009), Wang and Zhang (2012), Tang et al. (2005), D’Orazio and Succi (2004), Chen et 
al. (1996)). Boundary conditions form part of the iterative procedure at every iterative step, 
and therefore, these conditions should be defined in view of the theoretical boundary 
conditions so that the macroscopic boundary conditions are satisfied. During the streaming 
process, the velocities of the particles which start moving towards the wall are cut by the 
boundary. Therefore, the role of boundary conditions is to systematically define these 
distribution functions on the surface by an appropriate approximation to best fit the geometry 
and problem under investigation. 

Boundary conditions in LBM were originally adopted from lattice-Gas-Automata (LGA), 
usually referred to as Bounce-back scheme, which was used to achieve no-slip velocity 
conditions at the wall, making it ideal to problems involving solid or moving boundaries 
conditions, flow over obstacles, flow in complex geometries etc. The rules in bounce-back 
scheme are framed such a way that the distribution functions transmitted by a particle to a 
neighbouring node on the boundary while streaming are just returned back. It was established 
that this scheme has accuracy of only first order, whereas LBM has an accuracy of second 
order (Wang and Zhang (2012), Chen et al. (1996), He et al. (1997)). He et al. (1997) 
presented formula for slip-velocity in bounce-back rule. However, Suh et al. (2008) later 
revealed that the slip formula for the bounce-back algorithm presented in The et al. (1997) is 
erroneous and proposed a modified algorithm for better accuracy and stability of the LBM. 
They also showed that the numerical error of LBM occurs mainly due to slip-velocity at the 
solid boundaries. Many attempts have been made to improve accuracy of the bounce-back 
scheme. Inamuro et al. (1995) proposed boundary schemes which expressed the unknown 
proportion of particles on the boundary in terms of distribution function for wall and the 
counter slip velocity in the tangential direction. Maier et al. (1996) proposed to obtain the 
proportion of particles at the solid nodes using law of conservation of mass and momentum. 
Zou and He (1997) investigated pressure and velocity boundary conditions for 2-D and 3-D 
LBM-BGK models and proposed to define pressure and velocity on flow boundaries based on 
bounce-back of the non-equilibrium distribution function, which efficiently increased 
accuracy of the bounce-back scheme. Chen et al. (1996) proposed an extrapolation scheme to 
specify flow boundary conditions, in which they considered a layer inside the solid and 
distribution function was extrapolated from the boundary nodes and the first fluid node. Guo 
et al. (2002) proposed a first-order non-equilibrium extrapolation scheme in which the 
distribution function at the boundary node is specified in terms of the distribution function at 
the nodes in the solid wall according to the standard streaming-collision rule of the LBM. 

The Lattice Boltzmann Equation for flow in Porous Media 

The LBM-BGK equation for incompressible flow in porous media is given by 

/. (x + e j dt,t + dt ) - f t (x,t) = — — + d t F (1) 

T 

where z is the relaxation parameter, f) (x, / ) represents the proportion of particles at time t 
positioned at x moving with a velocity e, . The second term on the right hand side represents 
the effect of porous material to fluid flow which is included by adding total body force due to 
porous media (Ergun (1952), Peng (2003), He and Luo (1997), Seta et al. (2006)), along with 
the inertial force that includes the viscous force and external force giving 
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A suitable choice for F in discrete form is discussed in [Guo and Zhao (2002), Mehrizi et 
al. (2013), Peng et al. (2004)) given by 
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2. BOUNDARY CONDITIONS FOR LBM 

2.1. Boundary Conditions for Velocity Distribution Function 

Bounce-back scheme is a node based approach to map solid walls to the lattice structure, 
wherein lattice nodes are divided into a solid, fluid and boundary nodes. The boundary nodes 
are governed by a bounce-back rule for collision, in which a particle hitting the boundary is 
bounced back in the opposite direction and its velocity is reversed. Consequently, the 
boundary shifts inside the fluid by half of the lattice length (Comubert et al. [21]), leading to 
scheme being also called as half way partial bounce-back scheme. 

If density at the boundary is known, the unknown distribution function at the boundary 
can be expressed as a linear combination of known distribution functions. The bounce back 
rule of the non-equilibrium distribution functions is written as (Zou and He (1997)) 

y neq j neq (4) 

where a and p are in opposite direction. The unknown distribution functions can be 
defined in a similar way by splitting the distribution function as a sum of equilibrium and 
non-equilibrium parts (Mohamad (2011), Zou and He (1997)). 

The treatment is independent of the direction and results in a first order error in velocity. 
The scheme is very simple to implement in LBM and satisfies law of conservation of mass 
and momentum. However, there are issues with this scheme in satisfying no-slip conditions. 
Consequently, the bounce-back scheme is also expressed as (Suh et al. (2008)) 

fa(yF)= fp{yF + dt)=f p {y,t) + ^{fp{y,i)] (5) 

An alternative to this scheme is where a particle hitting the boundary is reflected back in 
such a way that only the normal velocity component changes sign, whereas the tangential 
velocity component is not changed. This scheme is referred to as specular reflection scheme 
(Cornubert et al. (1991)), which produces a full slip boundary unlike the half way bounce- 
back scheme. 

Another approach to predict distribution function at the boundary is by extrapolation and 
interpolation. This scheme involves extrapolating the distribution functions on the boundary 
in terms of the distribution functions on the first and second fluid nodes. For example, the 
unknown distribution functions on the west wall are given by 

fi U = l) = 2 ft {.x = 2) - f. (x = 3) (6) 

The unknown distribution functions on the other boundaries are defined in a similar way. 
Extrapolation schemes have been found to increase accuracy when the boundaries are curved 
or inclined. These schemes are found to increase the accuracy of the algorithm. One reason 
being that is these schemes carry out the collision step at the nodes on the solid wall. For 
example, the distribution function at the solid node, given by (Guo et al. (2002)), expresses it 


http://www.iaeme.com/IJMET/index.asp 



15 


editor@iaeme.com 



Shear-Thinning Fluid Flow in Porous Media: A Study of Boundary Behaviour Using Lattice 

Boltzmann Method 


as the sum of the equilibrium component and the non-equilibrium component with the 
equilibrium component of the distribution function defined using Eq. 3.6. 

fi (>’ = m) = f, ec ' (y = m) + f" eq {y = m) ( 7 ) 

The velocity component derived with the extrapolation scheme is given by (Suh et al. 
(2008)) 

u s = y l u 1 + y 2 u 2 (8) 


Where y and y 2 are constants that depend on the distance of the boundary wall to the 
nearest grid line towards the fluid side (w ) , given by 
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The non-equilibrium part is given by 


f’’ eq {x = n, y = m) = A f, neq (x = 1, y = l) + (l - A)f] neq {x = 2 ,y = 2) 

where the parameter Tis given by 
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Higher values of w enhances accuracy of the algorithm but results in less stability which 
can be decreased by decreasing the value of y . Besides extrapolation, interpolation can also 
be used to approximate distribution function on the boundary. However, this scheme needs a 
minimum number of nodes and is non-local in nature. A wide number of interpolation 
schemes are proposed in literature, the simplest among which can be found in (Yu et al. 
(2003)). These schemes are similar to bounce-back scheme, the only difference being that the 
value of distribution function depends on the node position on the wall. Yu et al. (2003) 
proposed the unknown distribution functions as 

f a (- x = Uy = u) = [(l - w)fg{x = i,y = 1 ,t) + wfg{x = i,y = l,t + dt )]+ ... 

i + w (13) 

+ •••- f a {x = i,y = 2,t) 

1 + w 


Both the schemes have second-order accuracy and the slip velocity increases with the 
relaxation parameter. In case of interpolation schemes, the relaxation parameter increases 
quadratically whereas for extrapolation the relaxation parameter increases linearly (Suh et al. 
(2008)). Extrapolation scheme, therefore, has an advantage over other schemes as the 
magnitude of error will be comparatively less for higher values of relaxation parameter. 

In case of repeating flow situations where the flow properties in upper portion of the 
geometry are identical to the lower portion, simulation can be simplified by ignoring the 
repeated flows to reduce computational time. Periodic boundary conditions can be used to 
model flows when a part of flow in upper portion of the geometry is identical to the part of 
flow in lower portion so that the boundary values on one side of the region are implemented 
by transferring the densities with velocity components from other side. These conditions are 
usually useful in problems involving bulk flow; flow in symmetric geometries and in case the 
surface effects are of very less interest. In case there is a perfect symmetry between flows 
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along a line, it is convenient to model the flow on only one side of the line, since the flow 
properties on one side is a reflection of the other side. 


2.2. Boundary Conditions for Thermal Distribution Function 

In general, boundary conditions for energy distribution function can be identical to that of 
velocity distribution function, the simplest approach being the bounce-back scheme and the 
second order extrapolation rule. The idea is to express heat flux in terms of kinetic moments 
of particle distribution function. Consequently, a number of ideas for implementing thermal 
boundary conditions have been proposed, some of which are briefed in subsequent sections. 
This approach is identical to the bounce-back rule of non-equilibrium parts for velocity 
distribution function. The unknown distribution function at the boundary can be expressed as 
a linear combination of known distribution functions. The bounce-back rule of the non- 
equilibrium distribution functions (Scheme I) is written as 

*r -‘i/r =-fcr -«w) (14> 


Where a and /? are in opposite direction, g t (x,t) is the energy distribution function and e 


is the internal energy given by e = — . 


Other approach presented by D’Orazio and Succi (2004) (Scheme II) suggested thermal 
energy density and heat flux as kinetic moments of separate energy distribution function 
g t (x,t) by doubling the degrees of freedom rather than supplying heat flux in terms of ki netic 
moments of velocity distribution function. For instance, the distribution functions on right 
wall were defined as 


8i = f (l - T p ) , g 5 = ^-r«(l + 3(« + V)) , g 8 = ^-T eq (1 + 3 (u - v)) 
y 3o 3o 

Where 

Tp = So + Si + g-i + §4 + <?6 + Si and T eq = -j — — ^ . (15) 

1 + 3 u 

Use of these boundary conditions is widely seen in simulating forced convection in porous 
media. However, difficulties in implementing these conditions were observed in case of 
natural convection. 

Thermal boundary conditions are also defined based on node weights in which the 
unknown distribution functions at the wall are specified at the boundaries as per Succi (2007) 
(Scheme III) 

8a= T w{ w a +w p)- 8p ( 16 ) 


These types of boundary conditions were observed to be very effective in simulating 
natural convection porous media. Though, a number of other ideas to implement thermal 
boundary conditions can be found in literature, it will depend on the application to identify the 
best suitable boundary conditions based on required fluid-structure interaction and geometry 
of the boundary. Description of some of these boundary conditions can be found in Succi 
(2007), Ho et al. (2009) and Chun and Ladd (2007). 
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3. INFLUENCE OF BOUNDARY CONDITIONS ON FLOW REGIMES 
IN POROUS MEDIA 

Boundary conditions discussed in Section 3 were applied to simulate fluid flow in porous 
media and the influence of boundary treatment was studied on flow regimes. This was done 
by simulating forced convection in porous media by a combination of two boundary 
treatments for velocity distribution function and three boundary treatments for energy 
distribution function (Scheme I-III). Consequently, influence of these boundary conditions 
was studied on velocity profiles and the rate of heat transfer. 




(a) (b) 

Figure 1 Comparison of velocity profiles at center of the cavity for Da=10 4 , e = 0.99 (a) as published 
in Guo and Zhao [18] (b) as computed in present investigation. 

3.1. Problem Description 

A square geometry containing porous medium was considered with top lid moving from left 
to right with a uniform velocity U = 0.1. The velocities at all other wall are assumed to be 
zero. A uniform fluid density /? = 1 .0 is imposed initially. Porosity e is set to 0.1 and the 
relaxation parameter r is assumed to be 0.8. Half-way bounce back boundary conditions are 
used in the numerical simulation. Numerical simulation is carried out for values of Da 
between 10 5 and 10 2 at Re = 10 and Re = 100. A MATLAB code was developed for a 
lOOxlOOmesh grid to solve the governing equations stated in Eq. 3.5-3.11 and Eq. 3.33. The 
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numerical technique was verified by comparing the results for Da=10 4 and e = 0.99 at 
Re = 400 and Re = 1000with Guo and Zhao [20]. The results are shown in Fig 1 and shows 
fine agreement with established results. 


4. RESULTS AND DISCUSSIONS 


This problem was to study the impact of boundary behavior on the velocity profiles and rate 
of heat transfer by considering forced convection influenced by a temperature difference 
between west and east wall of the cavity. Carreau-Yasuda model is used to represent non- 
Newtonian fluids given by Sochi (2010) 


M = Mo + iMo 



n — 1 
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1 + £t 
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V 

V J 
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where jU 0 and // x are viscosities at zero and infinite shear rate, respectively. The 
parameters a and b along with ju 0 and //^should be appropriately defined for numerical 

simulation to converge which depend on index n and material properties of the porous media. 
Simulation is carried out by using different types of boundary conditions each for velocity 
distribution function and energy distribution function. The boundary behaviour for velocity 
distribution function was considered based on half-way partial bounce back scheme and the 
bounce back scheme based on non-equilibrium distribution function, while behaviour for 
energy distribution function was considered based on scheme I, scheme II and scheme ITT. 
The momentum equation and energy equation are not coupled for forced convection; 
therefore, the energy equation can be solved after the momentum equations are solved to 
obtain velocity components. Both the schemes considered for velocity distribution function 
produced the same results as seen in Figure 1 . However, significant difference in rate of heat 
transfer was observed based on the scheme employed for energy distribution function. This 
was observed based on the variation in Nu depending on the thermal boundary conditions 
employed in the algorithm. Results were obtained for shear-thinning fluid by considering 
n = 0.5 at Da = l(r 2 and Re = 0.1 for various values of porosity, which are presented in Table 
1 . 

Table 1 Comparison of Nu values for two boundary schemes 



£ 

0.1 

0.4 

0.7 

Scheme I 

3.9717 

4.4537 

5.1441 

Scheme II 

18.0768 

19.0784 

21.0792 

Scheme III 

1.9802 

1.9889 

1.9902 


The difference in Nu is because of the fact that scheme II is able to include the viscous 
heating effects much better as compared to other schemes. This may be due to the ability of 
scheme express the heat flux in terms of only the first moment of energy distribution function. 
However, this scheme has stability issues when used in case of natural convection. Therefore, 
scheme I and III are a preferred choice for problems involving natural convection in porous 
media. In case of forced convection, scheme II produces a stronger heat transfer due to a 
better particle boundary interaction. The relative increase in the rate of heat transfer 
depending on porosity was observed in case of scheme I with the other two schemes 
producing a similar margin of increase between them. These observations are useful while 
simulating non-Newtonian forced and natural convection (non-isothermal flows) in porous 
media. 
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